## This code creates Figure 6 (a) and Figure 6 (b).

sample <- colonial_towns_shp[is.na(colonial_towns_shp$houses)==FALSE,]

sample$replaced_label <- ifelse(sample$replaced==0,"No","Yes") # X axis labels

# Figure 6 (a)

png(file=paste0("Plots/descriptive-all.png",sep=""), width=6, height=6, pointsize=9, units="in", res=600)

twowayplot <- ggplot(sample@data,
            aes(x=as.factor(replaced_label),fill=as.factor(coercive))) +
            geom_bar(width =.3,position = position_stack(reverse = T)) + 
            theme_bw() + scale_fill_grey() + xlab("New headman appointed before 1783") +
            ylab("Frequency") + labs(fill="Numer of colonial police stations") +
  theme(text = element_text(size = 16), axis.text.x = element_text(size=15),legend.position = "bottom")

twowayplot

dev.off()

# Simple t-test for the full sample

t.test(sample$coercive[sample$replaced==1],sample$coercive[sample$replaced==0])

# Figure 6 (b)

png(file=paste0("Plots/descriptive-500.png",sep=""), width=6, height=6, pointsize=9, units="in", res=600)

twowayplot_smallpop <- ggplot(sample@data[sample$houses<=500,],
                     aes(x=as.factor(replaced_label),fill=as.factor(coercive))) +
  geom_bar(width =.3,position = position_stack(reverse = T)) + 
  theme_bw() + scale_fill_grey() + xlab("New headman appointed before 1783") +
  ylab("Frequency") + labs(fill="Numer of colonial police stations") +
  theme(text = element_text(size = 16), axis.text.x = element_text(size=16),legend.position = "bottom")

twowayplot_smallpop

dev.off()


# Simple t-test for the observations with less than 500 houses only

t.test(sample$coercive[sample$replaced==1 & sample$houses<=500],sample$coercive[sample$replaced==0 & sample$houses<=500])
